
%% Either run these scripts or load their output files
%setup; 
load data_structures;
%Main_mpec
load est_point;
%counterfactualSolve
%load counterFact_out


%%
% Assume correct rhos are read in from est_point as rhohat
rho_model = e_rhohat;

rho_matrix_ger    = reshape(rho_model(1:(m.num_firms_ger+1)*m.num_pro_ger)',m.num_firms_ger+1,m.num_pro_ger)' ;
rho_matrix_dnk    = reshape(rho_model((m.num_firms_ger+1)*m.num_pro_ger+1:end)',m.num_firms_dnk+1,m.num_pro_dnk)' ;
win_dnk_in_ger_fit = sum(rho_matrix_ger(:,2:m.num_for_ger+1),2);
win_dnk_in_dnk_fit = sum(rho_matrix_dnk(:,2:end-1,:),2);

win_dnk_in_ger = sum(m.y_matrix_ger(:,2:m.num_for_ger+1),2);
win_dnk_in_dnk = sum(m.y_matrix_dnk(:,2:end-1,:),2);

rho_matrix_ger_nf    = reshape(e_rho_nF(1:(m_nFC.num_firms_ger+1)*m_nFC.num_pro_ger)',m_nFC.num_firms_ger+1,m_nFC.num_pro_ger)' ;
rho_matrix_dnk_nf    = reshape(e_rho_nF((m_nFC.num_firms_ger+1)*m_nFC.num_pro_ger+1:end)',m_nFC.num_firms_dnk+1,m_nFC.num_pro_dnk)' ;
win_dnk_in_ger_nf    = sum(rho_matrix_ger_nf(:,2:m_nFC.num_for_ger+1),2);
win_dnk_in_dnk_nf    = sum(rho_matrix_dnk_nf(:,2:end-m_nFC.num_for_dnk,:),2);

rho_matrix_ger_nb    = reshape(e_rho_nB(1:(m_nFC.num_firms_ger+1)*m_nFC.num_pro_ger)',m_nFC.num_firms_ger+1,m_nFC.num_pro_ger)' ;
rho_matrix_dnk_nb    = reshape(e_rho_nB((m_nFC.num_firms_ger+1)*m_nFC.num_pro_ger+1:end)',m_nFC.num_firms_dnk+1,m_nFC.num_pro_dnk)' ;
win_dnk_in_ger_nb    = sum(rho_matrix_ger_nb(:,2:m_nFC.num_for_ger+1),2);
win_dnk_in_dnk_nb    = sum(rho_matrix_dnk_nb(:,2:end-m_nFC.num_for_dnk,:),2);

x = [m.dist_bor_ger ; -m.dist_bor_dnk ];
y = [win_dnk_in_ger ; win_dnk_in_dnk ];
y_fit = [win_dnk_in_ger_fit; win_dnk_in_dnk_fit];
y_nf = [win_dnk_in_ger_nf; win_dnk_in_dnk_nf];
y_nb = [win_dnk_in_ger_nb; win_dnk_in_dnk_nb];


%Scale x to be between -1 and 1
x_s =  2*((x - min(x))./(max(x) - min(x)))-1;
pts = [-1:.005:1]';
ptsX = ((pts+1)./2).*(max(x) - min(x)) + min(x);

order = 5;
N = length(y);
X = zeros(N,order);
P = zeros(length(pts), order);
X(:,1) = ones(N,1);
P(:,1) = ones(length(pts),1);
X(:,2) = x_s;
P(:,2) = pts;
for c = 3:order
   X(:,c) = 2*x_s.*X(:,c-1) - X(:,c-2);
   P(:,c) = 2*pts.*P(:,c-1) - P(:,c-2);
end
%Here we calculate a 'simple' border effect'
X(:,order+1) = (x>=0);
P(:,order+1) = (ptsX >=0);

%Do ''full interaction''
%X = [X X.*(repmat((x>=0),1,order))];
%P = [P P.*(repmat((ptsX>=0), 1, order))];

b_hat      = (X'*X)\X'*y;
b_hat_fit  = (X'*X)\X'*y_fit;
b_hat_nf   = (X'*X)\X'*y_nf;
b_hat_nb   = (X'*X)\X'*y_nb;


f = P*b_hat;
f_fit = P*b_hat_fit;
f_nf = P*b_hat_nf;
f_nb = P*b_hat_nb;


plot(ptsX,f,'r');
hold on;
plot(ptsX,f_fit,'b');
plot(ptsX,f_nf,'g');
plot(ptsX,f_nb,'c');

saveas(gcf, 'border_fig', 'pdf');

